
library(tidyverse)
library(broom)
library(fixest)
library(ggthemes)
library(modelsummary)

## Main data
# source("code/01-importdata.R")


## Workflow
# 1. fixest::feols() models for main results
# 2. Tidying them into a data.frame
# 3. Making results table
# 4. Plotting
# 5. Interpretation

#### Main regression models ####

## Models with country-specific time trends

#### Run baseline models as from temperature shocks ####

### Annual count of disasters, for three lags
m_temp_annual1 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ annual_average_temp_lag1 + annual_W_lag1 |
                           iso3c[year] + year, # trends
                         data = analysis_df) %>% 
  as.list()

m_temp_annual2 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ annual_average_temp_lag2 + annual_W_lag2 |
                          iso3c[year] + year, # trends
                        data = analysis_df) %>% 
  as.list()

m_temp_annual3 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ annual_average_temp_lag3 + annual_W_lag3 |
                          iso3c[year] + year, # trends
                        data = analysis_df) %>% 
  as.list()

m_temp_season1 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ hottest_season_temp_lag1 + season_W_lag1 |
                          iso3c[year] + year, # trends
                        data = analysis_df) %>% 
  as.list()

m_temp_season2 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ hottest_season_temp_lag2 + season_W_lag2 |
                          iso3c[year] + year, # trends
                        data = analysis_df) %>% 
  as.list()

m_temp_season3 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ hottest_season_temp_lag3 + season_W_lag3 |
                          iso3c[year] + year, # trends
                        data = analysis_df) %>% 
  as.list()

temp_list <- do.call(c, list(m_temp_annual1, m_temp_annual2, m_temp_annual3,
                             m_temp_season1, m_temp_season2, m_temp_season3))

# models_trends <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
#                scaled_elec_renew, scaled_fit, scaled_quota,
#                scaled_delegates, scaled_gcg,
#                scaled_cf_total_provided, scaled_cf_principal_provided,
#                scaled_cf_total_received, scaled_cf_principal_received)
#              ~ sw(annual_average_temp_lag1, annual_average_temp_lag2,
#                   annual_average_temp_lag3,
#                   hottest_season_temp_lag1, hottest_season_temp_lag2,
#                   hottest_season_temp_lag3) |
#                iso3c[year] + year, # trends
#              data = analysis_df)
# list_trends <- as.list(models_trends)

results_trends <-  setNames(as.data.frame(matrix(NA, nrow = 66, ncol = 10)), 
                       c("outcome", "treatment", "beta", "std_error",
                         "t_stat", "p_value", "nobs", "n_countries",
                         "min_year", "max_year"))

for (i in 1:66){
  results_trends[i, 1] <- temp_list[[i]]$fml %>% gsub("\\~", "", .) %>% .[2] # outcome
  results_trends[i, 2] <- temp_list[[i]]$fml %>% gsub("\\~", "", .) %>% .[3] # treatment
  results_trends[i, 3] <- temp_list[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(2) # beta
  results_trends[i, 4] <- temp_list[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(3) # std error
  results_trends[i, 5] <- temp_list[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(4) # t-stat
  results_trends[i, 6] <- temp_list[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(5) # p-val
  results_trends[i, 7] <- temp_list[[i]]$nobs # number of observations
  results_trends[i, 8] <- length(unique(temp_list[[i]]$fixef_id$iso3c)) # number of countries
  results_trends[i, 9] <- min(attr(temp_list[[i]]$fixef_id$year, "fixef_names")) # start year of analysis
  results_trends[i, 10] <- max(attr(temp_list[[i]]$fixef_id$year, "fixef_names")) # end year of analysis
}
results_trends$treatment <- str_sub(results_trends$treatment, start = 1L, end = -17L)

# results_trends 


## Results table
results_trends %>% 
  filter(treatment == "annual_average_temp_lag1" |
           treatment == "hottest_season_temp_lag1") %>% 
  mutate(years = paste(min_year, max_year, sep="--")) %>% 
  select(outcome, treatment, beta, p_value, n_countries, years) %>% 
  mutate_at(vars(c(beta, p_value)), ~sprintf("%.3f", round(.,3))) %>% 
  pivot_wider(names_from = treatment, values_from = beta:p_value) %>% 
  mutate(labels = case_when(outcome=="scaled_climate_laws" ~ "Climate laws",
                            outcome=="scaled_fit" ~ "Feed-in tariff",
                            outcome=="scaled_quota" ~ "Renewables quota",
                            outcome=="scaled_elec_renew" ~ "Renewable electricity generation",
                            outcome=="scaled_ecp_zeros" ~ "Emissions-weighted carbon price",
                            outcome=="scaled_delegates" ~ "COP delegates",
                            outcome=="scaled_gcg" ~ "Climate institutional memberships",
                            outcome=="scaled_cf_total_provided" ~ "Climate finance provided",
                            outcome=="scaled_cf_principal_provided" ~ "Climate finance (principal) provided",
                            outcome=="scaled_cf_total_received" ~ "Climate finance received",
                            outcome=="scaled_cf_principal_received" ~ "Climate finance (principal) received")) %>% 
  select(labels, beta_annual_average_temp_lag1, p_value_annual_average_temp_lag1,
         beta_hottest_season_temp_lag1, p_value_hottest_season_temp_lag1,
         n_countries, years) %>% 
  knitr::kable(format = "html") %>% 
  kableExtra::save_kable(file = "text/table-global-oneyearlag-html.html")
  # knitr::kable(.,format = "latex", 
  #              linesep = "",
  #              booktabs = T,
  #              caption = "Average treatment effect of temperature shocks on climate policy")
# -> paste into manuscript

#### Plot results from trends models

# Create a labels object
labels_vector <- results_trends %>% 
  mutate(labels = case_when(outcome=="scaled_climate_laws" ~ "Climate laws",
                            outcome=="scaled_fit" ~ "Feed-in tariff",
                            outcome=="scaled_quota" ~ "Renewables quota",
                            outcome=="scaled_elec_renew" ~ "Renewable electricity\n generation",
                            outcome=="scaled_ecp_zeros" ~ "Emissions-weighted\n carbon price",
                            outcome=="scaled_delegates" ~ "COP delegates",
                            outcome=="scaled_gcg" ~ "Climate institutional\n memberships",
                            outcome=="scaled_cf_total_provided" ~ "Climate finance\n provided",
                            outcome=="scaled_cf_principal_provided" ~ "Climate finance\n (principal) provided",
                            outcome=="scaled_cf_total_received" ~ "Climate finance\n received",
                            outcome=="scaled_cf_principal_received" ~ "Climate finance\n (principal) received")) %>%
  select(labels)

## Plot temperature and lags in two facets
results_trends %>% 
  # Create lag ID
  mutate(Temperature = case_when(str_detect(treatment, 'lag1')==TRUE ~ "One-year lag",
                                 str_detect(treatment, 'lag2')==TRUE ~ "Two-year lag",
                                 str_detect(treatment, 'lag3')==TRUE ~ "Three-year lag")) %>% 
  # Create facet ID
  mutate(Facet = case_when(str_detect(treatment, 'season')==TRUE ~ "Hottest season",
                           str_detect(treatment, 'annual')==TRUE ~ "Annual")) %>% 
  # Create an index to organize results
  group_by(treatment) %>% 
  mutate(n = row_number()) %>% 
  ungroup() %>% 
  group_by(outcome, Facet) %>% 
  mutate(m = row_number()) %>% 
  ungroup() %>% 
  mutate(index = case_when(m == 1 ~ n - 0.2,
                           m == 2 ~ n - 0.0,
                           m == 3 ~ n + 0.2)) %>% 
  # Get confidence intervals
  mutate(ci_lower = beta - 1.96*std_error,
         ci_lower_correction = beta - 2.84*std_error,
         ci_upper = beta + 1.96*std_error,
         ci_upper_correction = beta + 2.84*std_error) %>%
  # Plot
  ggplot(., aes(beta, index, shape = Temperature)) +
  facet_grid(cols = vars(Facet)) +
  geom_segment(aes(x = ci_lower_correction, xend = ci_upper_correction,
                   y = index, yend = index),
               size = 1, color = "#79b321") + 
  # scale_color_manual(values = c("#79b321", "#578018", "#456613", "#344d0e")) +
  geom_segment(aes(x = ci_lower, xend = ci_upper, 
                   y = index, yend = index), 
               size = 2, color = "#456613") +
  geom_point(size = 3.5, aes(shape = Temperature)) +
  scale_shape_manual(breaks=c("One-year lag","Two-year lag","Three-year lag"),
                     values = c(16,3,4)) +
  scale_color_manual(values = c("black", "black", "black")) +
  geom_vline(xintercept = 0, linetype = 1, colour = "black") +
  scale_y_continuous(trans = "reverse",
                     breaks = c(1:11),
                     labels = unique(labels_vector$labels)) +
  labs(x="Regression coefficient and confidence intervals", 
       y="") +
  geom_vline(xintercept = 0) +
  theme(panel.spacing = unit(1, "lines")) +
  theme_tufte(base_family = "Helvetica") +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill=NA),
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.minor.x = element_line(colour = "grey90"),
        # panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=16))
ggsave("text/coefplot_lags_trends.pdf", units = "in", height = 8, width = 10)



modelsummary(m_temp_annual1, stars = T)
modelsummary(m_temp_season1, stars = T)
